
# clear environment
rm(list=ls())

# make sure working directory is set to folder within which this script is saved
getwd()

# load data
neg_bias <- read.csv("neg_bias_analysis_file.csv")

#### Combined Models ####

full_data <- subset(neg_bias, neg_bias$task != "loss") #remove loss aversion under ambiguity data

## Main
fit_polid <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)

fit_conservative <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)

fit_social <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)

fit_econ <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)

## Engagement interactions
fit_int_polid <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + as.factor(task), data = full_data)

fit_int_conservative <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + as.factor(task), data = full_data)

fit_int_social <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + as.factor(task), data = full_data)

fit_int_econ <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + as.factor(task), data = full_data)

## Personality traits
fit_open <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)

fit_nfc <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)

fit_conserve <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)

fit_auth <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + as.factor(task), data = full_data)


#### Loss Aversion under Ambiguity Models (Lucid 2018) ####

a_data <- subset(neg_bias, neg_bias$task == "loss")

## Main
fit_polid_a <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_conservative_a <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_social_a <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_econ_a <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

## Engagement interactions
fit_int_polid_a <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = a_data)

fit_int_conservative_a <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = a_data)

fit_int_social_a <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = a_data)

fit_int_econ_a <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = a_data)

## Personality traits
fit_open_a <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_nfc_a <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_conserve_a <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_auth_a <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

## Restricted measure (1/3 to 3)
fit_polid_3_a <- lm(right01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_conservative_3_a <- lm(conservative01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_social_3_a <- lm(social01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)

fit_econ_3_a <- lm(econ01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = a_data)


#### Lex Models #### 

l_data <- subset(neg_bias, neg_bias$task == "lex")

## Main
fit_polid_l <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

fit_conservative_l <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

fit_social_l <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

fit_econ_l <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

## Engagement interactions
fit_int_polid_l <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)

fit_int_conservative_l <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)

fit_int_social_l <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)

fit_int_econ_l <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)

## Personality traits
fit_open_l <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

fit_nfc_l <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

fit_conserve_l <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

fit_auth_l <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)

#### Lex Robust Models (only correct responses)

lex_tab <- array(NA, c(12,3))

## Main
fit_polid_l_cor <- lm(right01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)
lex_tab[1,] <- cbind(paste(format(round(coef(fit_polid_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_polid_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_polid_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_polid_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_polid_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_l_cor <- lm(conservative01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)
lex_tab[2,] <- cbind(paste(format(round(coef(fit_conservative_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_conservative_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_conservative_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_conservative_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_conservative_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_l_cor <- lm(social01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)
lex_tab[3,] <- cbind(paste(format(round(coef(fit_social_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_social_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_social_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_social_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_social_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_l_cor <- lm(econ01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)
lex_tab[4,] <- cbind(paste(format(round(coef(fit_econ_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_econ_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_econ_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_econ_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_econ_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

## Engagement interactions
fit_int_polid_l_cor <- lm(right01 ~ dscore_correct_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)
lex_tab[5,] <- cbind(paste(format(round(coef(fit_int_polid_l_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_polid_l_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_polid_l_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_polid_l_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_polid_l_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_l_cor <- lm(conservative01 ~ dscore_correct_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)
lex_tab[6,] <- cbind(paste(format(round(coef(fit_int_conservative_l_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_conservative_l_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_conservative_l_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_conservative_l_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_conservative_l_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_social_l_cor <- lm(social01 ~ dscore_correct_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)
lex_tab[7,] <- cbind(paste(format(round(coef(fit_int_social_l_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_social_l_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_social_l_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_social_l_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_social_l_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_econ_l_cor <- lm(econ01 ~ dscore_correct_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = l_data)
lex_tab[8,] <- cbind(paste(format(round(coef(fit_int_econ_l_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_econ_l_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_econ_l_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_econ_l_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_econ_l_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

## Personality traits
fit_open_l_cor <- lm(open01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)
lex_tab[9,] <- cbind(paste(format(round(coef(fit_open_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_open_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_open_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_open_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_open_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_l_cor <- lm(nclose01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)
lex_tab[10,] <- cbind(paste(format(round(coef(fit_nfc_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_nfc_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_nfc_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_nfc_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_nfc_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_l_cor <- lm(conserve01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 +
    income01 + unemp, data = l_data)
lex_tab[11,] <- cbind(paste(format(round(coef(fit_conserve_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conserve_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conserve_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conserve_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conserve_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_l_cor <- lm(auth01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = l_data)
lex_tab[12,] <- cbind(paste(format(round(coef(fit_auth_l_cor)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_auth_l_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_auth_l_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_auth_l_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_auth_l_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

rownames(lex_tab) <- c("Identity","All Issues","Social","Economic",
                        "Identity Interaction","All Issues Interaction","Social Interaction","Econ Interaction",
                        "(Low) Openness","Need for Closure","Conservation","Authoritarianism")
colnames(lex_tab) <- c("B (5%, 95%)","t", "p")
write.csv(lex_tab, "Tables/Supp K_Lex.csv")

#### Flank Models ####

f_data <- subset(neg_bias, neg_bias$task == "flank")

## Main
fit_polid_f <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

fit_conservative_f <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

fit_social_f <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

fit_econ_f <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

## Engagement interactions
fit_int_polid_f <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = f_data)

fit_int_conservative_f <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = f_data)

fit_int_social_f <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = f_data)

fit_int_econ_f <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = f_data)

## Personality traits
fit_open_f <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

fit_nfc_f <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

fit_conserve_f <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

fit_auth_f <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = f_data)

#### Flank Robust Models (only correct responses)

flk_tab <- array(NA, c(20,3))

## Main
fit_polid_f_cor <- lm(right01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
                        unemp, data = f_data)
flk_tab[1,] <- cbind(paste(format(round(coef(fit_polid_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_polid_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_polid_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_polid_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_polid_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_f_cor <- lm(conservative01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
                               unemp, data = f_data)
flk_tab[2,] <- cbind(paste(format(round(coef(fit_conservative_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_conservative_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_conservative_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_conservative_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_conservative_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_f_cor <- lm(social01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
                         unemp, data = f_data)
flk_tab[3,] <- cbind(paste(format(round(coef(fit_social_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_social_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_social_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_social_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_social_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_f_cor <- lm(econ01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = f_data)
flk_tab[4,] <- cbind(paste(format(round(coef(fit_econ_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_econ_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_econ_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_econ_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_econ_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

## Engagement interactions
fit_int_polid_f_cor <- lm(right01 ~ dscore_correct_std * engage01 + age01 + female + black +
                            hisp + educ01 + income01 + unemp, data = f_data)
flk_tab[5,] <- cbind(paste(format(round(coef(fit_int_polid_f_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_polid_f_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_polid_f_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_polid_f_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_polid_f_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_f_cor <- lm(conservative01 ~ dscore_correct_std * engage01 + age01 + female + black +
                                   hisp + educ01 + income01 + unemp, data = f_data)
flk_tab[6,] <- cbind(paste(format(round(coef(fit_int_conservative_f_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_conservative_f_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_conservative_f_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_conservative_f_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_conservative_f_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_social_f_cor <- lm(social01 ~ dscore_correct_std * engage01 + age01 + female + black +
                             hisp + educ01 + income01 + unemp, data = f_data)
flk_tab[7,] <- cbind(paste(format(round(coef(fit_int_social_f_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_social_f_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_social_f_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_social_f_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_social_f_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_econ_f_cor <- lm(econ01 ~ dscore_correct_std * engage01 + age01 + female + black +
                           hisp + educ01 + income01 + unemp, data = f_data)
flk_tab[8,] <- cbind(paste(format(round(coef(fit_int_econ_f_cor)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_econ_f_cor)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_econ_f_cor)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_econ_f_cor)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_econ_f_cor)$coef[11,4]), 2), nsmall = 3, trim = T))

## Personality traits
fit_open_f_cor <- lm(open01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = f_data)
flk_tab[9,] <- cbind(paste(format(round(coef(fit_open_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_open_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_open_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_open_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_open_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_f_cor <- lm(nclose01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
                      unemp, data = f_data)
flk_tab[10,] <- cbind(paste(format(round(coef(fit_nfc_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_nfc_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_nfc_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_nfc_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_nfc_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_f_cor <- lm(conserve01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 +
                           income01 + unemp, data = f_data)
flk_tab[11,] <- cbind(paste(format(round(coef(fit_conserve_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conserve_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conserve_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conserve_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conserve_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_f_cor <- lm(auth01 ~ dscore_correct_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = f_data)
flk_tab[12,] <- cbind(paste(format(round(coef(fit_auth_f_cor)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_auth_f_cor)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_auth_f_cor)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_auth_f_cor)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_auth_f_cor)$coef[2,4]), 2), nsmall = 3, trim = T))

#### Flank McLean et al. models
fit_polid_f_pos <- lm(right01 ~ flk_pos_only + age01 + female + black + hisp + educ01 + income01 +
                        unemp, data = f_data)
flk_tab[13,] <- cbind(paste(format(round(coef(fit_polid_f_pos)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_polid_f_pos)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_polid_f_pos)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_polid_f_pos)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_polid_f_pos)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_f_pos <- lm(conservative01 ~ flk_pos_only + age01 + female + black + hisp + educ01 + income01 +
                               unemp, data = f_data)
flk_tab[14,] <- cbind(paste(format(round(coef(fit_conservative_f_pos)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conservative_f_pos)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conservative_f_pos)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conservative_f_pos)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conservative_f_pos)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_f_pos <- lm(social01 ~ flk_pos_only + age01 + female + black + hisp + educ01 + income01 +
                         unemp, data = f_data)
flk_tab[15,] <- cbind(paste(format(round(coef(fit_social_f_pos)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_social_f_pos)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_social_f_pos)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_social_f_pos)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_social_f_pos)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_f_pos <- lm(econ01 ~ flk_pos_only + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = f_data)
flk_tab[16,] <- cbind(paste(format(round(coef(fit_econ_f_pos)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_econ_f_pos)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_econ_f_pos)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_econ_f_pos)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_econ_f_pos)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_polid_f_neg <- lm(right01 ~ flk_neg_only + age01 + female + black + hisp + educ01 + income01 +
                        unemp, data = f_data)
flk_tab[17,] <- cbind(paste(format(round(coef(fit_polid_f_neg)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_polid_f_neg)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_polid_f_neg)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_polid_f_neg)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_polid_f_neg)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_f_neg <- lm(conservative01 ~ flk_neg_only + age01 + female + black + hisp + educ01 + income01 +
                               unemp, data = f_data)
flk_tab[18,] <- cbind(paste(format(round(coef(fit_conservative_f_neg)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conservative_f_neg)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conservative_f_neg)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conservative_f_neg)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conservative_f_neg)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_f_neg <- lm(social01 ~ flk_neg_only + age01 + female + black + hisp + educ01 + income01 +
                         unemp, data = f_data)
flk_tab[19,] <- cbind(paste(format(round(coef(fit_social_f_neg)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_social_f_neg)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_social_f_neg)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_social_f_neg)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_social_f_neg)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_f_neg <- lm(econ01 ~ flk_neg_only + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = f_data)
flk_tab[20,] <- cbind(paste(format(round(coef(fit_econ_f_neg)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_econ_f_neg)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_econ_f_neg)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_econ_f_neg)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_econ_f_neg)$coef[2,4]), 2), nsmall = 3, trim = T))

rownames(flk_tab) <- c("Identity","All Issues","Social","Economic",
                       "Identity Interaction","All Issues Interaction","Social Interaction","Econ Interaction",
                       "(Low) Openness","Need for Closure","Conservation","Authoritarianism",
                       "Identity (Pos)","All Issues (Pos)","Social (Pos)","Econ (Pos)",
                       "Identity (Neg)","All Issues (Neg)","Social (Neg)","Econ (Neg)")
colnames(flk_tab) <- c("B (5%, 95%)","t", "p")
write.csv(flk_tab, "Tables/Supp K_Flank.csv")

#### Toubia Loss Aversion Models (Qualtrics Data) ####

t_data <- subset(neg_bias, neg_bias$task == "toubia")

## Main
fit_polid_t <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

fit_conservative_t <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

fit_social_t <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

fit_econ_t <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

## Engagement interactions
fit_int_polid_t <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = t_data)

fit_int_conservative_t <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = t_data)

fit_int_social_t <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = t_data)

fit_int_econ_t <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = t_data)

## Personality traits
fit_open_t <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

fit_nfc_t <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

fit_conserve_t <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

fit_auth_t <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = t_data)

### Robustness checks

toub_tab <- array(NA, c(24,3))

### Restricted measure (1/3 to 3)

## ideology
fit_polid_3_t <- lm(right01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                      unemp, data = t_data)
toub_tab[1,] <- cbind(paste(format(round(coef(fit_polid_3_t)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_polid_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_polid_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_polid_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_polid_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_3_t <- lm(conservative01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                             unemp, data = t_data)
toub_tab[2,] <- cbind(paste(format(round(coef(fit_conservative_3_t)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conservative_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conservative_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conservative_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conservative_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_3_t <- lm(social01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = t_data)
toub_tab[3,] <- cbind(paste(format(round(coef(fit_social_3_t)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_social_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_social_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_social_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_social_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_3_t <- lm(econ01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = t_data)
toub_tab[4,] <- cbind(paste(format(round(coef(fit_econ_3_t)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_econ_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_econ_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_econ_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_econ_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

## engagement interaction
fit_int_polid_3_t <- lm(right01 ~ kwloss_3_std * engage01 + age01 + female + black +
                          hisp + educ01 + income01 + unemp, data = t_data)
toub_tab[5,] <- cbind(paste(format(round(coef(fit_int_polid_3_t)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_polid_3_t)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_polid_3_t)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_polid_3_t)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_polid_3_t)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_3_t <- lm(conservative01 ~ kwloss_3_std * engage01 + age01 + female + black +
                                 hisp + educ01 + income01 + unemp, data = t_data)
toub_tab[6,] <- cbind(paste(format(round(coef(fit_int_conservative_3_t)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_conservative_3_t)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_conservative_3_t)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_conservative_3_t)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_conservative_3_t)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_social_3_t <- lm(social01 ~ kwloss_3_std * engage01 + age01 + female + black +
                           hisp + educ01 + income01 + unemp, data = t_data)
toub_tab[7,] <- cbind(paste(format(round(coef(fit_int_social_3_t)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_social_3_t)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_social_3_t)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_social_3_t)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_social_3_t)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_econ_3_t <- lm(econ01 ~ kwloss_3_std * engage01 + age01 + female + black +
                         hisp + educ01 + income01 + unemp, data = t_data)
toub_tab[8,] <- cbind(paste(format(round(coef(fit_int_econ_3_t)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_econ_3_t)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_econ_3_t)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_econ_3_t)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_econ_3_t)$coef[11,4]), 2), nsmall = 3, trim = T))

## personality
fit_open_3_t <- lm(open01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = t_data)
toub_tab[9,] <- cbind(paste(format(round(coef(fit_open_3_t)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_open_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_open_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_open_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_open_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_3_t <- lm(nclose01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                    unemp, data = t_data)
toub_tab[10,] <- cbind(paste(format(round(coef(fit_nfc_3_t)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_nfc_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_nfc_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_nfc_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_nfc_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_3_t <- lm(conserve01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                         unemp, data = t_data)
toub_tab[11,] <- cbind(paste(format(round(coef(fit_conserve_3_t)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_conserve_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_conserve_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_conserve_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_conserve_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_3_t <- lm(auth01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = t_data)
toub_tab[12,] <- cbind(paste(format(round(coef(fit_auth_3_t)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_auth_3_t)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_auth_3_t)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_auth_3_t)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_auth_3_t)$coef[2,4]), 2), nsmall = 3, trim = T))

### Include additional PT parameters as controls

## ideology
fit_polid_t_ro <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp + alpha + sigma, data = t_data)
toub_tab[13,] <- cbind(paste(format(round(coef(fit_polid_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_polid_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_polid_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_polid_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_polid_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_t_ro <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                              unemp + alpha + sigma, data = t_data)
toub_tab[14,] <- cbind(paste(format(round(coef(fit_conservative_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_conservative_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_conservative_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_conservative_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_conservative_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_t_ro <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                        unemp + alpha + sigma, data = t_data)
toub_tab[15,] <- cbind(paste(format(round(coef(fit_social_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_social_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_social_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_social_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_social_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_t_ro <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                      unemp + alpha + sigma, data = t_data)
toub_tab[16,] <- cbind(paste(format(round(coef(fit_econ_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_econ_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_econ_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_econ_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_econ_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

## engagement interaction
fit_int_polid_t_ro <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
                           hisp + educ01 + income01 + unemp + alpha + sigma, data = t_data)
toub_tab[17,] <- cbind(paste(format(round(coef(fit_int_polid_t_ro)[13],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_polid_t_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_polid_t_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_polid_t_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_polid_t_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_t_ro <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
                                  hisp + educ01 + income01 + unemp + alpha + sigma, data = t_data)
toub_tab[18,] <- cbind(paste(format(round(coef(fit_int_conservative_t_ro)[13],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_conservative_t_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_conservative_t_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_conservative_t_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_conservative_t_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

fit_int_social_t_ro <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
                            hisp + educ01 + income01 + unemp + alpha + sigma, data = t_data)
toub_tab[19,] <- cbind(paste(format(round(coef(fit_int_social_t_ro)[13],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_social_t_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_social_t_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_social_t_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_social_t_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

fit_int_econ_t_ro <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
                          hisp + educ01 + income01 + unemp + alpha + sigma, data = t_data)
toub_tab[20,] <- cbind(paste(format(round(coef(fit_int_econ_t_ro)[13],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_econ_t_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_econ_t_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_econ_t_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_econ_t_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

## personality
fit_open_t_ro <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                      unemp + alpha + sigma, data = t_data)
toub_tab[21,] <- cbind(paste(format(round(coef(fit_open_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_open_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_open_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_open_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_open_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_t_ro <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp + alpha + sigma, data = t_data)
toub_tab[22,] <- cbind(paste(format(round(coef(fit_nfc_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_nfc_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_nfc_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_nfc_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_nfc_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_t_ro <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp +
                          educ01 + income01 + unemp + alpha + sigma, data = t_data)
toub_tab[23,] <- cbind(paste(format(round(coef(fit_conserve_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_conserve_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_conserve_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_conserve_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_conserve_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_t_ro <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                      unemp + alpha + sigma, data = t_data)
toub_tab[24,] <- cbind(paste(format(round(coef(fit_auth_t_ro)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_auth_t_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_auth_t_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_auth_t_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_auth_t_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

rownames(toub_tab) <- rep(c("Identity","All Issues","Social","Economic",
                            "Identity Interaction","All Issues Interaction","Social Interaction","Econ Interaction",
                            "(Low) Openness","Need for Closure","Conservation","Authoritarianism"), 2)
colnames(toub_tab) <- c("B (5%, 95%)","t", "p")
write.csv(toub_tab, "Tables/Supp K_Toubia.csv")

#### Tanaka Loss Aversion Models (Lucid 2019 Data) ####

s_data <- subset(neg_bias, neg_bias$task == "summer19")

## Main
fit_polid_s <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

fit_conservative_s <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

fit_social_s <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

fit_econ_s <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

## Engagement interactions
fit_int_polid_s <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = s_data)

fit_int_conservative_s <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = s_data)

fit_int_social_s <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = s_data)

fit_int_econ_s <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp, data = s_data)

## Personality traits
fit_open_s <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

fit_nfc_s <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

fit_conserve_s <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

fit_auth_s <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp, data = s_data)

### Robustness checks

tan_tab <- array(NA, c(24,3))

### Restricted measure (1/3 to 3)

## ideology
fit_polid_3_s <- lm(right01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                      unemp, data = s_data)
tan_tab[1,] <- cbind(paste(format(round(coef(fit_polid_3_s)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_polid_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_polid_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_polid_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_polid_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_3_s <- lm(conservative01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                             unemp, data = s_data)
tan_tab[2,] <- cbind(paste(format(round(coef(fit_conservative_3_s)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_conservative_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_conservative_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_conservative_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_conservative_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_3_s <- lm(social01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = s_data)
tan_tab[3,] <- cbind(paste(format(round(coef(fit_social_3_s)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_social_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_social_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_social_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_social_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_3_s <- lm(econ01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = s_data)
tan_tab[4,] <- cbind(paste(format(round(coef(fit_econ_3_s)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_econ_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_econ_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_econ_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_econ_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

## engagement interaction
fit_int_polid_3_s <- lm(right01 ~ kwloss_3_std * engage01 + age01 + female + black +
                          hisp + educ01 + income01 + unemp, data = s_data)
tan_tab[5,] <- cbind(paste(format(round(coef(fit_int_polid_3_s)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_polid_3_s)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_polid_3_s)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_polid_3_s)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_polid_3_s)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_3_s <- lm(conservative01 ~ kwloss_3_std * engage01 + age01 + female + black +
                                 hisp + educ01 + income01 + unemp, data = s_data)
tan_tab[6,] <- cbind(paste(format(round(coef(fit_int_conservative_3_s)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_conservative_3_s)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_conservative_3_s)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_conservative_3_s)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_conservative_3_s)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_social_3_s <- lm(social01 ~ kwloss_3_std * engage01 + age01 + female + black +
                           hisp + educ01 + income01 + unemp, data = s_data)
tan_tab[7,] <- cbind(paste(format(round(coef(fit_int_social_3_s)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_social_3_s)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_social_3_s)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_social_3_s)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_social_3_s)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_econ_3_s <- lm(econ01 ~ kwloss_3_std * engage01 + age01 + female + black +
                         hisp + educ01 + income01 + unemp, data = s_data)
tan_tab[8,] <- cbind(paste(format(round(coef(fit_int_econ_3_s)[11],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_int_econ_3_s)[11,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_int_econ_3_s)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_int_econ_3_s)$coef[11,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_int_econ_3_s)$coef[11,4]), 2), nsmall = 3, trim = T))

## personality
fit_open_3_s <- lm(open01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = s_data)
tan_tab[9,] <- cbind(paste(format(round(coef(fit_open_3_s)[2],3), nsmall = 3, trim = T), " (", 
                           format(round(confint(fit_open_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                           format(round(confint(fit_open_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                     format(round(as.vector(summary(fit_open_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                     format(round(as.vector(summary(fit_open_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_3_s <- lm(nclose01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                    unemp, data = s_data)
tan_tab[10,] <- cbind(paste(format(round(coef(fit_nfc_3_s)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_nfc_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_nfc_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_nfc_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_nfc_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_3_s <- lm(conserve01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                         unemp, data = s_data)
tan_tab[11,] <- cbind(paste(format(round(coef(fit_conserve_3_s)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conserve_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conserve_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conserve_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conserve_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_3_s <- lm(auth01 ~ kwloss_3_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = s_data)
tan_tab[12,] <- cbind(paste(format(round(coef(fit_auth_3_s)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_auth_3_s)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_auth_3_s)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_auth_3_s)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_auth_3_s)$coef[2,4]), 2), nsmall = 3, trim = T))

### Include additional PT parameters as controls

## ideology
fit_polid_s_ro <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + alpha + sigma, data = s_data)
tan_tab[13,] <- cbind(paste(format(round(coef(fit_polid_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_polid_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_polid_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_polid_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_polid_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_s_ro <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + alpha + sigma, data = s_data)
tan_tab[14,] <- cbind(paste(format(round(coef(fit_conservative_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conservative_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conservative_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conservative_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conservative_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_s_ro <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + alpha + sigma, data = s_data)
tan_tab[15,] <- cbind(paste(format(round(coef(fit_social_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_social_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_social_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_social_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_social_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_s_ro <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + alpha + sigma, data = s_data)
tan_tab[16,] <- cbind(paste(format(round(coef(fit_econ_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_econ_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_econ_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_econ_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_econ_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

## engagement interaction
fit_int_polid_s_ro <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + alpha + sigma, data = s_data)
tan_tab[17,] <- cbind(paste(format(round(coef(fit_int_polid_s_ro)[13],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_polid_s_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_polid_s_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_polid_s_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_polid_s_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_s_ro <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + alpha + sigma, data = s_data)
tan_tab[18,] <- cbind(paste(format(round(coef(fit_int_conservative_s_ro)[13],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_conservative_s_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_conservative_s_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_conservative_s_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_conservative_s_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

fit_int_social_s_ro <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + alpha + sigma, data = s_data)
tan_tab[19,] <- cbind(paste(format(round(coef(fit_int_social_s_ro)[13],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_social_s_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_social_s_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_social_s_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_social_s_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

fit_int_econ_s_ro <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
    hisp + educ01 + income01 + unemp + alpha + sigma, data = s_data)
tan_tab[20,] <- cbind(paste(format(round(coef(fit_int_econ_s_ro)[13],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_econ_s_ro)[13,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_econ_s_ro)[13,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_econ_s_ro)$coef[13,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_econ_s_ro)$coef[13,4]), 2), nsmall = 3, trim = T))

## personality
fit_open_s_ro <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + alpha + sigma, data = s_data)
tan_tab[21,] <- cbind(paste(format(round(coef(fit_open_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_open_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_open_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_open_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_open_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_s_ro <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + alpha + sigma, data = s_data)
tan_tab[22,] <- cbind(paste(format(round(coef(fit_nfc_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_nfc_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_nfc_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_nfc_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_nfc_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_s_ro <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp +
    educ01 + income01 + unemp + alpha + sigma, data = s_data)
tan_tab[23,] <- cbind(paste(format(round(coef(fit_conserve_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conserve_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conserve_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conserve_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conserve_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_s_ro <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
    unemp + alpha + sigma, data = s_data)
tan_tab[24,] <- cbind(paste(format(round(coef(fit_auth_s_ro)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_auth_s_ro)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_auth_s_ro)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_auth_s_ro)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_auth_s_ro)$coef[2,4]), 2), nsmall = 3, trim = T))

rownames(tan_tab) <- rep(c("Identity","All Issues","Social","Economic",
                       "Identity Interaction","All Issues Interaction","Social Interaction","Econ Interaction",
                       "(Low) Openness","Need for Closure","Conservation","Authoritarianism"), 2)
colnames(tan_tab) <- c("B (5%, 95%)","t", "p")
write.csv(tan_tab, "Tables/Supp K_Tanaka.csv")

#### BeanFest (Lucid 2020 Data) ####

b_data <- subset(neg_bias, neg_bias$task == "summer20")

## Main
fit_polid_b <- lm(right01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                    unemp, data = b_data)

fit_conservative_b <- lm(conservative01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = b_data)

fit_social_b <- lm(social01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = b_data)

fit_econ_b <- lm(econ01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)

## Engagement interactions
fit_int_polid_b <- lm(right01 ~ neg_bias_std * engage01 + age01 + female + black +
                        hisp + educ01 + income01 + unemp, data = b_data)

fit_int_conservative_b <- lm(conservative01 ~ neg_bias_std * engage01 + age01 + female + black +
                           hisp + educ01 + income01 + unemp, data = b_data)

fit_int_social_b <- lm(social01 ~ neg_bias_std * engage01 + age01 + female + black +
                         hisp + educ01 + income01 + unemp, data = b_data)

fit_int_econ_b <- lm(econ01 ~ neg_bias_std * engage01 + age01 + female + black +
                       hisp + educ01 + income01 + unemp, data = b_data)

## Personality traits
fit_open_b <- lm(open01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)

fit_nfc_b <- lm(nclose01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                  unemp, data = b_data)

fit_conserve_b <- lm(conserve01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = b_data)

fit_auth_b <- lm(auth01 ~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)

#### Robustness checks

bean_tab <- array(NA, c(24,3))

### Avoidance Behavior Models

## ideology
fit_polid_b_ab <- lm(right01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                    unemp, data = b_data)
bean_tab[1,] <- cbind(paste(format(round(coef(fit_polid_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_polid_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_polid_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_polid_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_polid_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_b_ab <- lm(conservative01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = b_data)
bean_tab[2,] <- cbind(paste(format(round(coef(fit_conservative_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_conservative_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_conservative_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_conservative_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_conservative_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_b_ab <- lm(social01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = b_data)
bean_tab[3,] <- cbind(paste(format(round(coef(fit_social_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_social_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_social_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_social_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_social_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_b_ab <- lm(econ01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)
bean_tab[4,] <- cbind(paste(format(round(coef(fit_econ_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_econ_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_econ_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_econ_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_econ_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

## engagement interactions
fit_int_polid_b_ab <- lm(right01 ~ avoid_bev_std * engage01 + age01 + female + black +
                        hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[5,] <- cbind(paste(format(round(coef(fit_int_polid_b_ab)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_polid_b_ab)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_polid_b_ab)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_polid_b_ab)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_polid_b_ab)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_b_ab <- lm(conservative01 ~ avoid_bev_std * engage01 + age01 + female + black +
                           hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[6,] <- cbind(paste(format(round(coef(fit_int_conservative_b_ab)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_conservative_b_ab)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_conservative_b_ab)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_conservative_b_ab)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_conservative_b_ab)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_social_b_ab <- lm(social01 ~ avoid_bev_std * engage01 + age01 + female + black +
                         hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[7,] <- cbind(paste(format(round(coef(fit_int_social_b_ab)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_social_b_ab)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_social_b_ab)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_social_b_ab)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_social_b_ab)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_econ_b_ab <- lm(econ01 ~ avoid_bev_std * engage01 + age01 + female + black +
                       hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[8,] <- cbind(paste(format(round(coef(fit_int_econ_b_ab)[11],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_int_econ_b_ab)[11,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_int_econ_b_ab)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_int_econ_b_ab)$coef[11,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_int_econ_b_ab)$coef[11,4]), 2), nsmall = 3, trim = T))

## personality
fit_open_b_ab <- lm(open01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)
bean_tab[9,] <- cbind(paste(format(round(coef(fit_open_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                            format(round(confint(fit_open_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                            format(round(confint(fit_open_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                      format(round(as.vector(summary(fit_open_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                      format(round(as.vector(summary(fit_open_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_b_ab <- lm(nclose01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                  unemp, data = b_data)
bean_tab[10,] <- cbind(paste(format(round(coef(fit_nfc_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_nfc_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_nfc_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_nfc_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_nfc_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_b_ab <- lm(conserve01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = b_data)
bean_tab[11,] <- cbind(paste(format(round(coef(fit_conserve_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_conserve_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_conserve_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_conserve_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_conserve_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_b_ab <- lm(auth01 ~ avoid_bev_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)
bean_tab[12,] <- cbind(paste(format(round(coef(fit_auth_b_ab)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_auth_b_ab)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_auth_b_ab)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_auth_b_ab)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_auth_b_ab)$coef[2,4]), 2), nsmall = 3, trim = T))

### Valence Weighting Models

## ideology
fit_polid_b_vw <- lm(right01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                    unemp, data = b_data)
bean_tab[13,] <- cbind(paste(format(round(coef(fit_polid_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_polid_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_polid_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_polid_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_polid_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conservative_b_vw <- lm(conservative01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = b_data)
bean_tab[14,] <- cbind(paste(format(round(coef(fit_conservative_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_conservative_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_conservative_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_conservative_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_conservative_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_social_b_vw <- lm(social01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                     unemp, data = b_data)
bean_tab[15,] <- cbind(paste(format(round(coef(fit_social_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_social_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_social_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_social_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_social_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_econ_b_vw <- lm(econ01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)
bean_tab[16,] <- cbind(paste(format(round(coef(fit_econ_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_econ_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_econ_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_econ_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_econ_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

## engagement interactions
fit_int_polid_b_vw <- lm(right01 ~ val_wt_std * engage01 + age01 + female + black +
                        hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[17,] <- cbind(paste(format(round(coef(fit_int_polid_b_vw)[11],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_polid_b_vw)[11,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_polid_b_vw)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_polid_b_vw)$coef[11,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_polid_b_vw)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_conservative_b_vw <- lm(conservative01 ~ val_wt_std * engage01 + age01 + female + black +
                           hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[18,] <- cbind(paste(format(round(coef(fit_int_conservative_b_vw)[11],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_conservative_b_vw)[11,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_conservative_b_vw)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_conservative_b_vw)$coef[11,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_conservative_b_vw)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_social_b_vw <- lm(social01 ~ val_wt_std * engage01 + age01 + female + black +
                         hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[19,] <- cbind(paste(format(round(coef(fit_int_social_b_vw)[11],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_social_b_vw)[11,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_social_b_vw)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_social_b_vw)$coef[11,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_social_b_vw)$coef[11,4]), 2), nsmall = 3, trim = T))

fit_int_econ_b_vw <- lm(econ01 ~ val_wt_std * engage01 + age01 + female + black +
                       hisp + educ01 + income01 + unemp, data = b_data)
bean_tab[20,] <- cbind(paste(format(round(coef(fit_int_econ_b_vw)[11],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_int_econ_b_vw)[11,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_int_econ_b_vw)[11,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_int_econ_b_vw)$coef[11,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_int_econ_b_vw)$coef[11,4]), 2), nsmall = 3, trim = T))

## personality
fit_open_b_vw <- lm(open01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)
bean_tab[21,] <- cbind(paste(format(round(coef(fit_open_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_open_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_open_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_open_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_open_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_nfc_b_vw <- lm(nclose01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                  unemp, data = b_data)
bean_tab[22,] <- cbind(paste(format(round(coef(fit_nfc_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_nfc_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_nfc_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_nfc_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_nfc_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_conserve_b_vw <- lm(conserve01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                       unemp, data = b_data)
bean_tab[23,] <- cbind(paste(format(round(coef(fit_conserve_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_conserve_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_conserve_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_conserve_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_conserve_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

fit_auth_b_vw <- lm(auth01 ~ val_wt_std + age01 + female + black + hisp + educ01 + income01 +
                   unemp, data = b_data)
bean_tab[24,] <- cbind(paste(format(round(coef(fit_auth_b_vw)[2],3), nsmall = 3, trim = T), " (", 
                             format(round(confint(fit_auth_b_vw)[2,1],3), nsmall = 3, trim = T), ", ", 
                             format(round(confint(fit_auth_b_vw)[2,2],3), nsmall = 3, trim = T), ")", sep=""), 
                       format(round(as.vector(summary(fit_auth_b_vw)$coef[2,3]), 2), nsmall = 3, trim = T), 
                       format(round(as.vector(summary(fit_auth_b_vw)$coef[2,4]), 2), nsmall = 3, trim = T))

rownames(bean_tab) <- rep(c("Identity","All Issues","Social","Economic",
                       "Identity Interaction","All Issues Interaction","Social Interaction","Econ Interaction",
                       "(Low) Openness","Need for Closure","Conservation","Authoritarianism"), 2)
colnames(bean_tab) <- c("B (5%, 95%)","t", "p")
write.csv(bean_tab, "Tables/Supp K_Bean.csv")

#### Figures ####

plot_info_builder <- function(model1, model2, model3, model4, model5, model6, coef_pos) {
  plot_info <- list()
  
  plot_info <- matrix(NA, nrow = 6, ncol = 3)
  plot_info[ ,1] <- c(coef(model1)[coef_pos], coef(model2)[coef_pos],
                      coef(model3)[coef_pos], coef(model4)[coef_pos],
                      coef(model5)[coef_pos], coef(model6)[coef_pos])
  
  plot_info[ ,2] <- c(confint(model1)[coef_pos,1], confint(model2)[coef_pos,1],
                      confint(model3)[coef_pos,1], confint(model4)[coef_pos,1],
                      confint(model5)[coef_pos,1], confint(model6)[coef_pos,1])
  
  plot_info[ ,3] <- c(confint(model1)[coef_pos,2], confint(model2)[coef_pos,2],
                      confint(model3)[coef_pos,2], confint(model4)[coef_pos,2],
                      confint(model5)[coef_pos,2], confint(model6)[coef_pos,2])
  
  return(plot_info)
}

## Political DVs Figure

plot_polid <- plot_info_builder(fit_polid, fit_polid_l, fit_polid_f, fit_polid_t, fit_polid_s, fit_polid_b, 2)
plot_conservative <- plot_info_builder(fit_conservative, fit_conservative_l, fit_conservative_f, fit_conservative_t, fit_conservative_s, fit_conservative_b, 2)
plot_social <- plot_info_builder(fit_social, fit_social_l,fit_social_f, fit_social_t, fit_social_s, fit_social_b, 2)
plot_econ <- plot_info_builder(fit_econ, fit_econ_l,fit_econ_f,fit_econ_t, fit_econ_s, fit_econ_b, 2)


pdf(file="Figures/Figure 1.pdf", height = 3.46, width = 3.46, pointsize=6, family = "Helvetica")

par(mfrow = c(2,2), mar=c(3,4,2,1))

plot(1:6, plot_polid[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_polid[1:6,2], 1:6, plot_polid[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Right-wing political identity")
box()

plot(1:6, plot_conservative[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_conservative[1:6,2], 1:6, plot_conservative[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="General conservatism")
box()

plot(1:6, plot_social[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_social[1:6,2], 1:6, plot_social[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Social conservatism")
box()

plot(1:6, plot_econ[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_econ[1:6,2], 1:6, plot_econ[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Economic conservatism")
box()

dev.off()


## Personality DVs Figure

plot_open <- plot_info_builder(fit_open, fit_open_l, fit_open_f, fit_open_t, fit_open_s, fit_open_b, 2)
plot_nfc <- plot_info_builder(fit_nfc, fit_nfc_l, fit_nfc_f, fit_nfc_t, fit_nfc_s, fit_nfc_b, 2)
plot_conserve <- plot_info_builder(fit_conserve, fit_conserve_l, fit_conserve_f, fit_conserve_t, fit_conserve_s, fit_conserve_b, 2)
plot_auth <- plot_info_builder(fit_auth, fit_auth_l, fit_auth_f, fit_auth_t, fit_auth_s, fit_auth_b, 2)


pdf(file="Figures/Figure 2.pdf", height = 3.46, width = 3.46, pointsize=6, family = "Helvetica")

par(mfrow = c(2,2), mar=c(3,4,2,1))

plot(1:6, plot_open[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_open[1:6,2], 1:6, plot_open[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="(Low) Openness to experience")
box()

plot(1:6, plot_nfc[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_nfc[1:6,2], 1:6, plot_nfc[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Need for closure")
box()

plot(1:6, plot_conserve[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_conserve[1:6,2], 1:6, plot_conserve[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Conservation v. openness")
box()

plot(1:6, plot_auth[1:6,1], ylab = "Marginal effect of negativity bias", xlab = "", axes = F, pch=16,
     ylim = c(-0.1,0.1), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.1,-0.05,0,0.05,0.1))
segments(1:6, plot_auth[1:6,2], 1:6, plot_auth[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Authoritarian childrearing values")
box()

dev.off()


## Engagement interaction figure

plot_int_polid <- plot_info_builder(fit_int_polid, fit_int_polid_l, fit_int_polid_f, fit_int_polid_t, fit_int_polid_s, fit_int_polid_b, 11)
plot_int_conservative <- plot_info_builder(fit_int_conservative, fit_int_conservative_l, fit_int_conservative_f, fit_int_conservative_t, fit_int_conservative_s, fit_int_conservative_b, 11)
plot_int_social <- plot_info_builder(fit_int_social, fit_int_social_l, fit_int_social_f, fit_int_social_t, fit_int_social_s, fit_int_social_b, 11)
plot_int_econ <- plot_info_builder(fit_int_econ, fit_int_econ_l,  fit_int_econ_f, fit_int_econ_t, fit_int_econ_s, fit_int_econ_b, 11)

# fix combined model
plot_int_polid[1, ] <- c(coef(fit_int_polid)[15], confint(fit_int_polid)[15,1],
                         confint(fit_int_polid)[15,2])
plot_int_conservative[1, ] <- c(coef(fit_int_conservative)[15], confint(fit_int_conservative)[15,1],
                          confint(fit_int_conservative)[15,2])
plot_int_social[1, ] <- c(coef(fit_int_social)[15], confint(fit_int_social)[15,1],
                         confint(fit_int_social)[15,2])
plot_int_econ[1, ] <- c(coef(fit_int_econ)[15], confint(fit_int_econ)[15,1],
                         confint(fit_int_econ)[15,2])


pdf(file="Figures/Figure 3.pdf", height = 3.46, width = 3.46, pointsize=6, family = "Helvetica")

par(mfrow = c(2,2), mar=c(3,4,2,1))

plot(1:6, plot_int_polid[1:6,1], ylab = "Coefficient for engagement interaction", xlab = "", axes = F, pch=16,
     ylim = c(-0.2,0.2), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.2,-0.1,0,0.1,0.2))
segments(1:6, plot_int_polid[1:6,2], 1:6, plot_int_polid[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Right-wing political identity")
box()

plot(1:6, plot_int_conservative[1:6,1], ylab = "Coefficient for engagement interaction", xlab = "", axes = F, pch=16,
     ylim = c(-0.2,0.2), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.2,-0.1,0,0.1,0.2))
segments(1:6, plot_int_conservative[1:6,2], 1:6, plot_int_conservative[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="General conservatism")
box()

plot(1:6, plot_int_social[1:6,1], ylab = "Coefficient for engagement interaction", xlab = "", axes = F, pch=16,
     ylim = c(-0.2,0.2), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.2,-0.1,0,0.1,0.2))
segments(1:6, plot_int_social[1:6,2], 1:6, plot_int_social[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Social conservatism")
box()

plot(1:6, plot_int_econ[1:6,1], ylab = "Coefficient for engagement interaction", xlab = "", axes = F, pch=16,
     ylim = c(-0.2,0.2), xlim = c(.5,6.5))
axis(1, at=c(1:6), labels = c("C","L","F","To","Ta","B"))
axis(2, at=c(-0.2,-0.1,0,0.1,0.2))
segments(1:6, plot_int_econ[1:6,2], 1:6, plot_int_econ[1:6,3], lty = 1, lwd = 0.5)
abline(h=0, lty = 2, lwd=0.5)
title(main="Economic conservatism")
box()

dev.off()

#### Table Maker ####

# Function to create tables from model objects
Table_Maker <- function(models, tab_names) {
  
  names1 <- c("Intercept", "Neg. Bias","Age (0-1)", "Female", "Black", "Hispanic", "Educ (0-1)", "Income (0-1)","Unemployed", "Task = Lex", "Task = Loss '19", "Task = Bean", "Task = Loss '14","N","Adj. R2")
  names2 <- c("Intercept", "Neg. Bias", "Pol. Engage (0-1)","Age (0-1)", "Female", "Black", "Hispanic", "Educ (0-1)", "Income (0-1)","Unemployed", "Task = Lex", "Task = Loss '19", "Task = Bean", "Task = Loss '14","Interaction","N","Adj. R2")
  names3 <- c("Intercept", "Neg. Bias","Age (0-1)", "Female", "Black", "Hispanic", "Educ (0-1)", "Income (0-1)","Unemployed","N","Adj. R2")
  names4 <- c("Intercept", "Neg. Bias", "Pol. Engage (0-1)","Age (0-1)", "Female", "Black", "Hispanic", "Educ (0-1)", "Income (0-1)","Unemployed","Interaction","N","Adj. R2")
  
  for (i in 1:7){
    
    x <- list()
    
    for (j in 1:4){
      
      B <- format(round(as.vector(summary(models[[j + (i-1)*4]])$coef[,1]), 2), nsmall = 2, trim = T)
      low <- format(round(as.vector(confint(models[[j + (i-1)*4]])[,1]), 2), nsmall = 2, trim = T)
      high <- format(round(as.vector(confint(models[[j + (i-1)*4]])[,2]), 2), nsmall = 2, trim = T)
      BCI <- paste(B, " (", low, ", ", high,")", sep = "")
      
      t <- format(round(as.vector(summary(models[[j + (i-1)*4]])$coef[,3]), 2), nsmall = 2, trim = T)
      p <- format(round(as.vector(summary(models[[j + (i-1)*4]])$coef[,4]), 2), nsmall = 2, trim = T)
      tp <- paste(t, " (", p, ")", sep="")

      
      n <- nobs(models[[j + (i-1)*4]])
      adj_r2 <- round(summary(models[[j + (i-1)*4]])$adj.r.squared, 2)
      x[[j]] <- cbind(append(BCI,c(n,adj_r2)), append(tp,c("","")))

    }
    
    tab <- do.call(cbind.data.frame, x)
    colnames(tab) <- rep(c("B (5%, 95%)","t (p)"), 4)
    
    if (length(B) == 13){
      rownames(tab) <- names1
    }else if (length(B) == 15){
      rownames(tab) <- names2
    }else if (length(B) == 9){
      rownames(tab) <- names3
    }else{
      rownames(tab) <- names4
    }

    write.csv(tab, paste("Tables/", tab_names[[i]], ".csv",  sep=""))

  }
}

# Figure 1 models and labels
fig1_models <- list(fit_polid, fit_conservative, fit_econ, fit_social, # combined models
                    fit_polid_l, fit_conservative_l, fit_econ_l, fit_social_l, # lex models 
                    fit_polid_f, fit_conservative_f, fit_econ_f, fit_social_f, # flanker models
                    fit_polid_t, fit_conservative_t, fit_econ_t, fit_social_t, # Toubia loss aversion models
                    fit_polid_s, fit_conservative_s, fit_econ_s, fit_social_s, # Tanaka loss aversion models
                    fit_polid_b, fit_conservative_b, fit_econ_b, fit_social_b, # beanfest models
                    fit_polid_a, fit_conservative_a, fit_econ_a, fit_social_a  # Loss under ambiguity models
)
fig1_names <- list("Ideo_Combined", # combined models
                   "Ideo_Lex", # lex models 
                   "Ideo_Flank", # flanker models
                   "Ideo_Toubia", # Toubia loss aversion models
                   "Ideo_Tanaka", # Tanaka loss aversion models
                   "Ideo_Bean", # beanfest models
                   "Ideo_Ambig" # Loss under ambiguity
)

# Figure 2 models and labels
fig2_models <- list(fit_open, fit_nfc, fit_conserve, fit_auth, # combined models
                    fit_open_l, fit_nfc_l, fit_conserve_l, fit_auth_l, # lex models 
                    fit_open_f, fit_nfc_f, fit_conserve_f, fit_auth_f, # flanker models
                    fit_open_t, fit_nfc_t, fit_conserve_t, fit_auth_t, # Toubia loss aversion models
                    fit_open_s, fit_nfc_s, fit_conserve_s, fit_auth_s, # Tanaka loss aversion models
                    fit_open_b, fit_nfc_b, fit_conserve_b, fit_auth_b, # beanfest models
                    fit_open_a, fit_nfc_a, fit_conserve_a, fit_auth_a # loss ambiguity
)
fig2_names <- list("Pers_Combined", # combined models
                   "Pers_Lex", # lex models 
                   "Pers_Flank", # flanker models
                   "Pers_Toubia", # Toubia loss aversion models
                   "Pers_Tanaka", # Tanaka loss aversion models
                   "Pers_Bean", # beanfest models
                   "Pers_Ambig" # Loss under ambiguity
)

# Figure 3 models and labels
fig3_models <- list(fit_int_polid, fit_int_conservative, fit_int_econ, fit_int_social, # combined models
                    fit_int_polid_l, fit_int_conservative_l, fit_int_econ_l, fit_int_social_l, # lex models 
                    fit_int_polid_f, fit_int_conservative_f, fit_int_econ_f, fit_int_social_f, # flanker models
                    fit_int_polid_t, fit_int_conservative_t, fit_int_econ_t, fit_int_social_t, # Toubia loss aversion models
                    fit_int_polid_s, fit_int_conservative_s, fit_int_econ_s, fit_int_social_s, # Tanaka loss aversion models
                    fit_int_polid_b, fit_int_conservative_b, fit_int_econ_b, fit_int_social_b, # beanfest models
                    fit_int_polid_a, fit_int_conservative_a, fit_int_econ_a, fit_int_social_a # loss ambiguity
)
fig3_names <- list("Int_Combined", # combined models
                   "Int_Lex", # lex models 
                   "Int_Flank", # flanker models
                   "Int_Toubia", # Toubia loss aversion models
                   "Int_Tanaka", # Tanaka loss aversion models
                   "Int_Bean", # beanfest models
                   "Int_Ambig" # Loss under ambiguity
)


## Make tables (saves all individual tables to 'Tables' sub-folder)

# Figure 1 tables
Table_Maker(fig1_models, fig1_names)

# Figure 2 tables
Table_Maker(fig2_models, fig2_names)

# Figure 3 tables
Table_Maker(fig3_models, fig3_names)


#### Correlation table of personality and ideology constructs ####

cors <- round(cor(data.frame(neg_bias$right01, neg_bias$social01, neg_bias$econ01,
                             neg_bias$conservative01, neg_bias$open01, neg_bias$nclose01,
                             neg_bias$conserve01, neg_bias$auth01), use = "complete.obs"), 2)
rownames(cors) <- colnames(cors) <- c("Identity","Social","Economic","General","Low Openness","NFC","Conservation","Auth")
write.csv(cors ,"Tables/Bivariate Correlations.csv")

#### Benjamini et al. Significance Adjustments ####

p <- c( # pull p-values
  
  summary(fit_polid)$coefficients[2,4], 
  summary(fit_conservative)$coefficients[2,4], 
  summary(fit_social)$coefficients[2,4], 
  summary(fit_econ)$coefficients[2,4], 
  summary(fit_nfc)$coefficients[2,4],
  summary(fit_open)$coefficients[2,4],
  summary(fit_conserve)$coefficients[2,4], 
  summary(fit_auth)$coefficients[2,4],
  
  summary(fit_polid_l)$coefficients[2,4], 
  summary(fit_conservative_l)$coefficients[2,4], 
  summary(fit_social_l)$coefficients[2,4], 
  summary(fit_econ_l)$coefficients[2,4], 
  summary(fit_nfc_l)$coefficients[2,4],
  summary(fit_open_l)$coefficients[2,4],
  summary(fit_conserve_l)$coefficients[2,4], 
  summary(fit_auth_l)$coefficients[2,4],
  
  summary(fit_polid_f)$coefficients[2,4], 
  summary(fit_conservative_f)$coefficients[2,4], 
  summary(fit_social_f)$coefficients[2,4], 
  summary(fit_econ_f)$coefficients[2,4], 
  summary(fit_nfc_f)$coefficients[2,4],
  summary(fit_open_f)$coefficients[2,4],
  summary(fit_conserve_f)$coefficients[2,4], 
  summary(fit_auth_f)$coefficients[2,4],
  
  summary(fit_polid_t)$coefficients[2,4], 
  summary(fit_conservative_t)$coefficients[2,4], 
  summary(fit_social_t)$coefficients[2,4], 
  summary(fit_econ_t)$coefficients[2,4], 
  summary(fit_nfc_t)$coefficients[2,4],
  summary(fit_open_t)$coefficients[2,4],
  summary(fit_conserve_t)$coefficients[2,4], 
  summary(fit_auth_t)$coefficients[2,4],
  
  summary(fit_polid_s)$coefficients[2,4], 
  summary(fit_conservative_s)$coefficients[2,4], 
  summary(fit_social_s)$coefficients[2,4], 
  summary(fit_econ_s)$coefficients[2,4], 
  summary(fit_nfc_s)$coefficients[2,4],
  summary(fit_open_s)$coefficients[2,4],
  summary(fit_conserve_s)$coefficients[2,4], 
  summary(fit_auth_s)$coefficients[2,4],
  
  summary(fit_polid_b)$coefficients[2,4], 
  summary(fit_conservative_b)$coefficients[2,4], 
  summary(fit_social_b)$coefficients[2,4], 
  summary(fit_econ_b)$coefficients[2,4], 
  summary(fit_nfc_b)$coefficients[2,4],
  summary(fit_open_b)$coefficients[2,4],
  summary(fit_conserve_b)$coefficients[2,4], 
  summary(fit_auth_b)$coefficients[2,4],
  
  summary(fit_int_polid)$coefficients[15,4], 
  summary(fit_int_conservative)$coefficients[15,4], 
  summary(fit_int_social)$coefficients[15,4], 
  summary(fit_int_econ)$coefficients[15,4], 
  
  summary(fit_int_polid_l)$coefficients[11,4], 
  summary(fit_int_conservative_l)$coefficients[11,4], 
  summary(fit_int_social_l)$coefficients[11,4], 
  summary(fit_int_econ_l)$coefficients[11,4], 
  
  summary(fit_int_polid_f)$coefficients[11,4], 
  summary(fit_int_conservative_f)$coefficients[11,4], 
  summary(fit_int_social_f)$coefficients[11,4], 
  summary(fit_int_econ_f)$coefficients[11,4], 
  
  summary(fit_int_polid_t)$coefficients[11,4], 
  summary(fit_int_conservative_t)$coefficients[11,4], 
  summary(fit_int_social_t)$coefficients[11,4], 
  summary(fit_int_econ_t)$coefficients[11,4], 
  
  summary(fit_int_polid_s)$coefficients[11,4], 
  summary(fit_int_conservative_s)$coefficients[11,4], 
  summary(fit_int_social_s)$coefficients[11,4], 
  summary(fit_int_econ_s)$coefficients[11,4], 
  
  summary(fit_int_polid_b)$coefficients[11,4], 
  summary(fit_int_conservative_b)$coefficients[11,4], 
  summary(fit_int_social_b)$coefficients[11,4], 
  summary(fit_int_econ_b)$coefficients[11,4]
  )

# row labels for table
col1 <- c( rep("combined",8), rep("lex",8),rep("flank",8), rep("toubia",8), rep("tanaka",8), rep("bean",8),
           rep("combined",4), rep("lex",4),rep("flank",4), rep("toubia",4), rep("tanaka",4), rep("bean",4))
col2 <- c(rep(c("Identity","General","Social","Econ","NFC","Low Openness","Conservation","Auth"), 6),
          rep(c("Identity","General","Social","Econ"), 6))

p_adj_1 <- p.adjust(p, method = 'BH', n = length(p))
rank(p)
tab1 <- cbind(col1, col2, rank(p), round(p,2), round(p_adj_1,2))
colnames(tab1) <- c("Dataset","Model","p-value rank","Original p-value","B-H adjusted p-value")
write.csv(tab1, "Tables/BH adjustment.csv")

p_adj_2 <- p.adjust(p, method = 'BY', n = length(p))
rank(p)
tab2 <- cbind(col1, col2, rank(p), round(p,2), round(p_adj_2,2))
colnames(tab2) <- c("Dataset","Model","p-value rank","Original p-value","B-Y adjusted p-value")
write.csv(tab2, "Tables/BY adjustment.csv")

#### Bayes Factors ####

# required packages
library(rstanarm)
library(bayestestR)
library(effectsize)
library(logspline)

# dependent variable list
dvs <- c('right01', 'conservative01', 'social01', 'econ01', 'nclose01', 'open01', 'conserve01', 'auth01')

# dataset list
datasets <- list(full_data, l_data, f_data, t_data, s_data, b_data)

# initialize vectors to hold Bayes factors and qualitative interpretations
bf <- qual <- c()

# model formulas
form1 <- "~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 + unemp + as.factor(task)"
form2 <- "~ neg_bias_std * engage01 + age01 + female + black + hisp + educ01 + income01 + unemp + as.factor(task)"
form3 <- "~ neg_bias_std + age01 + female + black + hisp + educ01 + income01 + unemp"
form4 <- "~ neg_bias_std * engage01 + age01 + female + black + hisp + educ01 + income01 + unemp"

# get estimates and Bayes factors
set.seed(8)
for(j in 1:length(datasets)) {
  
  for(i in 1:length(dvs)) {
    
    ## combined data
    if(j == 1) {
      
    # no interaction
      bayesfit1 <- stan_glm(
        formula = as.formula(paste(dvs[i], form1, sep="")), 
        data = datasets[[j]],
        prior = normal(0, 0.5, autoscale = T), 
        prior_intercept = normal(0.5, 0.5, autoscale = FALSE))
      BF_1 <- bayesfactor_parameters(bayesfit1, null = 0, direction = "right")
      bf <- append(bf, 1/exp(BF_1$log_BF[2]))
      qual <- append(qual, 
                     interpret_bf(1/exp(BF_1$log_BF[2]),
                                  include_value = TRUE, 
                                  rules = "raftery1995"))
      
    # interaction (only ideology DVs)
      if (i <= 4){
        bayesfit2 <- stan_glm(
          formula = as.formula(paste(dvs[i], form2, sep="")), 
          data = datasets[[j]],
          prior = normal(0, 0.5, autoscale = T),
          prior_intercept = normal(0.5, 0.5, autoscale = FALSE))
        BF_2 <- bayesfactor_parameters(bayesfit2, null = 0, direction = "right")
        bf <- append(bf, 1/exp(BF_2$log_BF[15]))
        qual <- append(qual, 
                       interpret_bf(1/exp(BF_2$log_BF[15]),
                                    include_value = TRUE, 
                                    rules = "raftery1995"))
      }
    }
    
    ## single studies
    else{
      
    # no interaction
      bayesfit1 <- stan_glm(
        formula = as.formula(paste(dvs[i], form3, sep="")), 
        data = datasets[[j]],
        prior = normal(0, 0.5, autoscale = T),
        prior_intercept = normal(0.5, 0.5, autoscale = FALSE))
      BF_1 <- bayesfactor_parameters(bayesfit1, null = 0, direction = "right")
      bf <- append(bf, 1/exp(BF_1$log_BF[2]))
      qual <- append(qual, 
                     interpret_bf(1/exp(BF_1$log_BF[2]),
                                  include_value = TRUE, 
                                  rules = "raftery1995"))
      
      # interaction (only ideology DVs)
      if (i <= 4){
        bayesfit2 <- stan_glm(
          formula = as.formula(paste(dvs[i], form4, sep="")), 
          data = datasets[[j]],
          prior = normal(0, 0.5, autoscale = T),
          prior_intercept = normal(0.5, 0.5, autoscale = FALSE))
        BF_2 <- bayesfactor_parameters(bayesfit2, null = 0, direction = "right")
        bf <- append(bf, 1/exp(BF_2$log_BF[11]))
        qual <- append(qual, 
                       interpret_bf(1/exp(BF_2$log_BF[11]),
                                    include_value = TRUE, 
                                    rules = "raftery1995"))
      }
    }
  }
}

# rownames for output
col1 <- c( rep("Combined",12), rep("Lexical",12),rep("Flanker",12), rep("Loss '14",12), rep("Loss '19",12), rep("BeanFest",12) )
col2 <- c(rep(c("Identity","Identity - Int",
                "General","General - Int",
                "Social","Social - Int",
                "Econ","Econ - Int",
                "NFC","Low Openness","Conservation","Auth"
                ) 
              , 6))

# write to csv
tab <- cbind(col1,col2,qual)
colnames(tab) <- c("Dataset","Model","BF and qualitative interpretation")
write.csv(tab, "Tables/Bayes Factors.csv")




